%% SCRIPT TO PROCESS CFRADIAL DOW DATA FOR MARSHALL FIRE TO PRODUCE 3D GRIDDED VOLUMES AND SURFACE MAPS
%produced by Neil P. Lareau summer 2022
%for questions contact: nlareau@unr.ed

close all
clear all


addpath('~/Dropbox/MATLAB/')
cd('~/Dropbox/MARSHALL_FIRE/')
re=6.37*10^6;%earth radius

%% Start by reading JSON file with fire perimtere data
fileName='marshall.json';
str = fileread(fileName); % dedicated for reading files as text
test=jsondecode(str);
figure(1);clf;
for aa=1:size(test.features.geometry.coordinates)
    figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot(pnow(:,1),pnow(:,2),'r')
end

%% List the DOW netcdf files this can be for deployment1 or deployment2, and you will need to change the target directory accordingly
files=dir('./deployment1/cfrad*.nc');

%% For the first pass through the data we need to group files into "volumes" defined as upward increments in the tilt angle of the radar
for ff=1:length(files);
    fname=strcat(files(ff).folder,'/',files(ff).name); %construct the filepath/filename
    azimuth=double(ncread(fname,'azimuth')); %parse out the azimuth data
    elevation=double(ncread(fname,'elevation')); %parse out the elevation data
    evlist(ff)=mode(elevation); %store the modal elevation angle for this sweep into a list (evlist)
end

%% Volumes are defined as having incrementally INCREASING elevations and at least 4 elevation steps
figure(1);clf;
nvals=1:ff; %total number of elevations in the list
dev=[0 diff(evlist)]; %this is the finite difference between elevation increments
sp1=subplot(2,1,1); cla; hold on;
plot(1:ff,evlist) %make a plot showing the ramp patterns in the elevation angle
hold on;
scatter(nvals(dev<-1),evlist(dev<-1),'*k') %these are down stepping periods, which mark the break points in the volume scanning
sp2=subplot(2,1,2);cla; hold on;
plot(nvals,dev)
hold on;
scatter(nvals(dev<-1),dev(dev<-1),'*k')

%% now find the dev<1 points as break points and group data into segement indices for the volumes
bidx=find(dev<-1); %indices for the break points
col=jet(50)
for bbb=1:(length(bidx)-1) 
    if bbb==1
        segment(bbb).indices=1:bidx(bbb);
    else
        segment(bbb).indices=bidx(bbb):bidx(bbb+1)
    end
    plot(sp1,segment(bbb).indices,evlist(segment(bbb).indices)) %this should plot each segement in a unique color and provides a visual check that this approach is working
end

%% Now loop over the number of segements, and add all the data from the corresponding files to one volume

% set up x,y,z 3d grid to house the data
xv=(-2*10^4):100:(4*10^4);
yv=(-1.6*10^4):100:(.5*10^4);
zv=[0:100:5000];
[xxx,yyy,zzz]=meshgrid(xv,yv,zv);
%initialize a master matrix for the relfectivity and time data
dbzmaster=nan(length(segment),length(yv),length(xv),length(zv));
timemaster=nan(length(segment),1);

%loop over the segements (i.e., volume scans)
for vv=1:length(segment)
    disp(vv./length(vv)) %this shows the fractional progress (e.g., .75 is 75% through the files)
    fileidx=segment(vv).indices; %the file #s in this segement
    allx=[]; %blank x points
    ally=[]; %blank y points
    allz=[]; %blank z points
    alldbz=[]; %blank dbz points

    for ff=1:length(fileidx) %loop over the relevant files
        fname=strcat(files(fileidx(ff)).folder,'/',files(fileidx(ff)).name); %construct the filename
        ncdisp(fname)%display the file name
        lats=ncread(fname,'latitude'); %parse out latitude
        lons=ncread(fname,'longitude'); %parse out longitude
        alt=ncread(fname,'altitude'); %parse out altitude
        vel=double(ncread(fname,'VEL')); %radial velocity 
        velF=double(ncread(fname,'VEL_F')); %unfolded radial velocoty
        dbzhc=double(ncread(fname,'DBZHC')); %horizontal polarization radar reflectivity
        azimuth=double(ncread(fname,'azimuth'));  %antena azimuth
        elevation=double(ncread(fname,'elevation')); %antenna elevation
        ranger=double(ncread(fname,'range')); %range gate distances
        stime=[ncread(fname,'time_coverage_start')']; %start of scan interval
        yyyymmdd=stime([1:4 6 7 9 10]); %pull out year, month,day
        HHMMSS=stime([12 13 15 16 18 19]); %pull out hour, minuteb, second
        timemaster(vv)=datenum(strcat(yyyymmdd,HHMMSS),'yyyymmddHHMMSS'); %construct matlab serial time

        clear xdist ydist dbznow

        for ev=1%:size(elevation) -> note this is hard coded to 1, but in some other similar scripts loops over elevations within a single file
            evnow=mean(elevation); % pull out the elevation for this scan.
            if evnow>0 %provided the elevaiton is positive
                ae=(4/3).*(6.37*10^6); %earth radius
                z=(((ranger.^2) + (ae.^2) + (2*ranger.*ae.*(sind(evnow)))).^(1/2))-ae; %height of the radar beam using 4/3 earth model
                hdist=  ae.*(asin((ranger.*cosd(evnow))./(ae+z)));% horizontal distance along earth's surface
                [azmt,hdstmt]=meshgrid(azimuth,hdist); % create a mesh of azimuths and along ground distances
                [~,zdist]=meshgrid(azimuth,z); %z grid
                xdist=sind(azmt).*hdstmt; %x component of the horizontal distance
                ydist=cosd(azmt).*hdstmt; %y comonent of the horizontal distance
                scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);%scale factor for longitude-to-x conversions
                xdistlon=xdist./(scalefactor); %conversion
                ydistlat=ydist./(111180); %conversion


                %% Produce Plan-view maps of the reflectivity and radial velocity data for each scan
                addpath('~/Dropbox/MATLAB/SGP_LIDAR_CODES/')

                figure(1);clf;set(gcf,'pos',[100 100 1200 500],'color','w')
                sp1=subplot(1,2,1);cla;hold on;
                colormap(sp1,rbmapper(30,-10)) %produces a divergent red blue colormap
                pcolor(xdist,ydist,vel);shading flat;caxis([-10 30]) %display the velocity data
                cbh=colorbar; %set colorbar on
                title(num2str(evnow),'fontsize',15) %title
                xlim([-1 3.5]*10^4); %xlimits to the plot
                ylim([-2 1]*10^4); % y-limits to the plot
                daspect([1 1 1]) %aspect ratio (Set to 1 here)
                sp2=subplot(1,2,2);cla; hold on;
                dbzhc=medfilt2(dbzhc); %smooth the reflectivity data with 3x3 median filter
                dbzhc(xdist<-3300)=nan; %remove values at large negative ranges (these are returns off the terrain to the west of Boulder)
                dbzhc(dbzhc<-5)=nan; %remove low values
                dbzhc(ydist>(0.2*10^4))=nan;%remove values far to the north
                dbzhc=medfilt2(dbzhc); %filter

                %% add padded edges around reflectivity features. This is very useful for the 3D renderings later in that it prevents routines from connecting reflectivity features across space that should not connect
                dbzhcx=imdilate(dbzhc,ones(10,10));
                dbzhc(~isinf(dbzhcx) & isnan(dbzhc))=-30;

              
                % display the buffered/padded reflecitivity data
                pcolor(xdist,ydist,(dbzhc));shading flat;caxis([-10 30])
                cbh=colorbar;
                colormap(sp2,jet(48))
                daspect([1 1 1])
                title(num2str(evnow),'fontsize',15)
                xlim([-1 3.5]*10^4);
                ylim([-2 1]*10^4);
                pause(.3)

                %% now add all the x,y,z,dbz data to point clouds to be used in the 3D interpolation
                allx=[allx(:); xdist(~isnan(dbzhc))];
                ally=[ally(:); ydist(~isnan(dbzhc))];
                allz=[allz(:); zdist(~isnan(dbzhc))];
                alldbz=[alldbz(:); dbzhc(~isnan(dbzhc))];

            end
        end
    end


    %% visualize scattered 3d data
    figure(188);clf;hold on;
    scatter3(allx(alldbz>10),ally(alldbz>10),allz(alldbz>10),10,alldbz(alldbz>10),'filled');
    caxis([0 30]);

    %% Perform the 3d interpolation from the scattered data to gridded data
    tic
    F1=TriScatteredInterp(allx(alldbz>-50),ally(alldbz>-50),allz(alldbz>-50),alldbz(alldbz>-50));
    toc

    tic
    dbzint=F1(xxx(:),yyy(:),zzz(:));
    toc
    dbzint=reshape(dbzint,size(xxx)); %need to reform as a 3D matrix
    dbzint(isnan(dbzint))=-50;%set nan values to large negatives for smooth surface rendering later
    
    dbzmaster(vv,:,:,:)=dbzint; %store interpolated data for output in dbzmaster

end


%% Output data
stationlat=nanmean(lats); %save the mean lat
stationlon=nanmean(lons); %save the mean long
stationalt=nanmean(alt); %save the mean elevation
%use this for deployment 1
save('marhsall_plume_d1.mat','dbzmaster','timemaster','xxx','yyy','zzz','stationlat','stationlon','stationalt','-v7.3')

%use this for deployment 2
% save('marhsall_plume_d2.mat','dbzmaster','timemaster','xxx','yyy','zzz','stationlat','stationlon','stationalt','-v7.3')


